WBE for SARS-CoV-2 in Switzerland: model development (A2)
Julien Riou
2023-06-08
source("setup.R")
ww1 = readRDS(fs::path("../",controls$savepoint,"ww1.rds"))
shapes = readRDS(fs::path("../",controls$savepoint,"shapes.rds"))We continue with model A5, now extending to multiple ARAs together. We start with NUTS-2 regions.
Mittelland
We start with the Mittelland region, that includes 24 ARAs.
# select NUTS2 region
select_NUTS2 = "Mittelland"
ww_reg = ww1 %>%
# log
dplyr::mutate(logvl=log(vl)) %>%
# select one ARA per NUTS-2
dplyr::filter(NUTS2_name==select_NUTS2) %>%
# replace zero values by 1
dplyr::mutate(vl=if_else(vl==0, 1, vl)) %>%
dplyr::filter(!is.na(vl)) %>%
# create indexes for INLA
dplyr::mutate(day1=day,
ara1=as.numeric(as.factor(ara_n)),
ara2=ara1)
# correspondence table
corr_reg = ww_reg %>%
group_by(ara_n,ara_id,ara_name,kt,pop,lab,lab_n,ara1,ara2) %>%
count() %>%
ungroup()
mw_100_desc_table(ww_reg) %>%
dplyr::mutate(across(everything(),as.character)) %>%
tidyr::gather() %>%
flextable::flextable(cwidth=c(4,4)) key | value |
|---|---|
Number of ARAs | 24 |
Number of laboratories | 3 |
Number of laboratory methods | 4 |
Number of measurements | 4554 |
Measurements below LOQ | 51 |
Measurements below LOD | 77 |
First | 2022-02-07 |
Last | 2023-05-10 |
Median viral concentration [gc/L] | 6e+04 (range: 0 to 4e+06) |
Median flow [m3/day] | 2e+04 (range: 4e+02 to 8e+05) |
Median viral load [gc/day/100,000] | 3e+12 (range: 1 to 1e+14) |
ara_name | Number of ARAs | Number of laboratories | Number of laboratory methods | Number of measurements | Measurements below LOQ | Measurements below LOD |
|---|---|---|---|---|---|---|
Aarwangen (Zala) | 1 | 1 | 1 | 137 | 0 | 2 |
Biel | 1 | 1 | 1 | 135 | 0 | 0 |
Burgdorf | 1 | 1 | 1 | 54 | 0 | 0 |
Colombier (La Saunerie) | 1 | 1 | 1 | 137 | 4 | 0 |
Delemont (Sede) | 1 | 1 | 1 | 279 | 0 | 0 |
Fribourg (Aele) | 1 | 1 | 1 | 329 | 0 | 8 |
Grenchen (Buerenamt) | 1 | 1 | 1 | 119 | 0 | 0 |
Grindelwald | 1 | 1 | 1 | 200 | 0 | 3 |
Interlaken | 1 | 1 | 1 | 134 | 0 | 0 |
La-Chaux-de-Fonds | 1 | 1 | 1 | 136 | 5 | 0 |
Laupen | 1 | 1 | 1 | 447 | 35 | 0 |
Lauterbrunnen | 1 | 1 | 1 | 143 | 0 | 1 |
Lyss | 1 | 1 | 1 | 133 | 0 | 2 |
Mittl. Emmental | 1 | 1 | 1 | 134 | 0 | 0 |
Moossee-Urtenenbach | 1 | 1 | 1 | 136 | 0 | 20 |
Neuchatel | 1 | 1 | 1 | 321 | 7 | 0 |
Porrentruy (Sepe) | 1 | 1 | 1 | 228 | 0 | 2 |
Region Bern | 1 | 1 | 1 | 280 | 0 | 11 |
Saanen | 1 | 1 | 1 | 151 | 0 | 25 |
Thunersee | 1 | 1 | 1 | 336 | 0 | 0 |
Vuippens (Ais) | 1 | 1 | 1 | 134 | 0 | 0 |
Winznau (Zv Olten) | 1 | 1 | 1 | 134 | 0 | 0 |
Worblental | 1 | 1 | 1 | 126 | 0 | 1 |
Zuchwil (Soloth.-Emme) | 1 | 1 | 1 | 191 | 0 | 2 |
ggplot(ww_reg) + geom_line(aes(x=date,y=logvl,colour=ara_name),alpha=.6) + scale_colour_discrete(guide="none") + labs(title="Viral load measurements")
Model A5.1: fixed effect
We add a geographical structure to model A5, starting by simply adding a fixed effect by ARA, assuming that the time trends of viral load across ARAs are identical. In Mittelland there is no change of methods, and there are 3 laboratories including 2 that only process one ARA each, so we leave these aspects out for now.
if(controls$rerun_models) {
ma5.1 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw1", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
ara_name,
data = ww_reg,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.1,file=paste0("../",controls$savepoint,"ma5.1.rds"))
} else {
ma5.1 = readRDS(file=paste0("../",controls$savepoint,"ma5.1.rds"))
}
summary(ma5.1)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 0.71, Running = 4.95, Post = 0.186, Total = 5.84
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 13.796 2319.820 -4654.827 13.797 4682.410 13.799 0
## ara_nameBiel 0.517 0.083 0.354 0.517 0.679 0.517 0
## ara_nameBurgdorf 0.023 0.111 -0.195 0.023 0.241 0.023 0
## ara_nameColombier (La Saunerie) 0.368 0.083 0.206 0.368 0.531 0.368 0
## ara_nameDelemont (Sede) 0.082 0.073 -0.061 0.082 0.225 0.082 0
## ara_nameFribourg (Aele) 0.120 0.070 -0.018 0.120 0.257 0.120 0
## ara_nameGrenchen (Buerenamt) 0.300 0.086 0.131 0.300 0.470 0.300 0
## ara_nameGrindelwald 0.821 0.078 0.669 0.821 0.973 0.821 0
## ara_nameInterlaken 0.787 0.084 0.623 0.787 0.952 0.787 0
## ara_nameLa-Chaux-de-Fonds -0.305 0.083 -0.468 -0.305 -0.142 -0.305 0
## ara_nameLaupen 0.698 0.069 0.564 0.698 0.833 0.698 0
## ara_nameLauterbrunnen 1.448 0.083 1.286 1.448 1.610 1.448 0
## ara_nameLyss 0.381 0.084 0.216 0.381 0.546 0.381 0
## ara_nameMittl. Emmental -0.259 0.084 -0.423 -0.259 -0.095 -0.259 0
## ara_nameMoossee-Urtenenbach 0.125 0.084 -0.040 0.125 0.290 0.125 0
## ara_nameNeuchatel 0.677 0.071 0.538 0.677 0.816 0.677 0
## ara_namePorrentruy (Sepe) 0.062 0.075 -0.085 0.062 0.209 0.062 0
## ara_nameRegion Bern -0.423 0.072 -0.565 -0.423 -0.282 -0.423 0
## ara_nameSaanen 0.257 0.082 0.097 0.257 0.417 0.257 0
## ara_nameThunersee 0.186 0.070 0.049 0.186 0.323 0.186 0
## ara_nameVuippens (Ais) 0.032 0.084 -0.132 0.032 0.196 0.032 0
## ara_nameWinznau (Zv Olten) 0.138 0.084 -0.026 0.138 0.302 0.138 0
## ara_nameWorblental 0.024 0.084 -0.141 0.024 0.189 0.024 0
## ara_nameZuchwil (Soloth.-Emme) 0.337 0.078 0.184 0.337 0.490 0.337 0
## weekend 0.007 0.030 -0.052 0.007 0.067 0.007 0
## hol 0.075 0.078 -0.078 0.075 0.228 0.075 0
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.18 0.044 2.09 2.178 2.266 2.18
## Precision for below_loq 9.14 5.036 2.42 8.126 21.708 6.06
## Precision for below_lod 0.00 0.000 0.00 0.000 0.000 0.00
## Precision for day 0.33 0.045 0.25 0.326 0.428 0.32
##
## Watanabe-Akaike information criterion (WAIC) ...: 267148.50
## Effective number of parameters .................: 542.74
##
## Marginal log-Likelihood: -134582.16
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')## exp(beta) 0.025quant 0.975quant
## weekend 1.01 0.95 1.07
## hol 1.08 0.93 1.26
## [1] "Relative viral load by ARA compared to reference 'Aarwangen (Zala)':"## exp(beta) 0.025quant 0.975quant
## ara_nameLauterbrunnen 4.25 3.62 5.00
## ara_nameGrindelwald 2.27 1.95 2.65
## ara_nameInterlaken 2.20 1.86 2.59
## ara_nameLaupen 2.01 1.76 2.30
## ara_nameNeuchatel 1.97 1.71 2.26
## ara_nameBiel 1.68 1.43 1.97
## ara_nameLyss 1.46 1.24 1.73
## ara_nameColombier (La Saunerie) 1.45 1.23 1.70
## ara_nameZuchwil (Soloth.-Emme) 1.40 1.20 1.63
## ara_nameGrenchen (Buerenamt) 1.35 1.14 1.60
## ara_nameSaanen 1.29 1.10 1.52
## ara_nameThunersee 1.20 1.05 1.38
## ara_nameWinznau (Zv Olten) 1.15 0.97 1.35
## ara_nameFribourg (Aele) 1.13 0.98 1.29
## ara_nameMoossee-Urtenenbach 1.13 0.96 1.34
## ara_nameDelemont (Sede) 1.09 0.94 1.25
## ara_namePorrentruy (Sepe) 1.06 0.92 1.23
## ara_nameVuippens (Ais) 1.03 0.88 1.22
## ara_nameBurgdorf 1.02 0.82 1.27
## ara_nameWorblental 1.02 0.87 1.21
## ara_nameMittl. Emmental 0.77 0.65 0.91
## ara_nameLa-Chaux-de-Fonds 0.74 0.63 0.87
## ara_nameRegion Bern 0.65 0.57 0.75We observe a large heterogeneity across ARAs, with highest relative viral load compared to reference in Lauterbrunnen, Interlaken and Grindelwald and lower relative viral load in Bern, La Chaux-de-Fonds and Emmental (note that viral load already accounts for population size).
Model A5.2: random intercept
We consider using a random intercept by ARA.
if(controls$rerun_models) {
ma5.2 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="iid"),
data = ww_reg,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.2,file=paste0("../",controls$savepoint,"ma5.2.rds"))
} else {
ma5.2 = readRDS(file=paste0("../",controls$savepoint,"ma5.2.rds"))
}
summary(ma5.2)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 1.98, Running = 69.5, Post = 0.145, Total = 71.6
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 21.259 328736.466 -649479.349 21.005 649523.338 20.517 0
## weekend -0.002 0.023 -0.047 -0.002 0.042 -0.002 0
## hol 0.103 0.066 -0.026 0.103 0.231 0.103 0
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.04 0.011 2.02 2.04 2.06 2.04
## Precision for below_loq 2.36 0.017 2.32 2.36 2.39 2.36
## Precision for below_lod 0.00 0.000 0.00 0.00 0.00 0.00
## Precision for day 0.01 0.000 0.01 0.01 0.01 0.01
## Precision for ara1 2.48 0.020 2.44 2.48 2.52 2.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 267513.60
## Effective number of parameters .................: 557.74
##
## Marginal log-Likelihood: -136432.44
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')## exp(beta) 0.025quant 0.975quant
## weekend 1.00 0.95 1.04
## hol 1.11 0.97 1.26
## [1] "Relative viral load by ARA compared to average:"ma5.2$summary.random$ara1 %>%
as_tibble() %>%
dplyr::transmute(ara1=ID,
`exp(beta)`=round(exp(mean),2),
`0.025quant`=round(exp(`0.025quant`),2),
`0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>%
dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>%
dplyr::arrange(-`exp(beta)`) %>%
select(-ara1) %>%
column_to_rownames(var="ara_name")## exp(beta) 0.025quant 0.975quant
## Lauterbrunnen 3.20 2.42 4.23
## Grindelwald 1.75 1.33 2.30
## Interlaken 1.67 1.26 2.20
## Laupen 1.50 1.15 1.95
## Neuchatel 1.49 1.14 1.94
## Biel 1.26 0.96 1.67
## Lyss 1.14 0.86 1.50
## Colombier (La Saunerie) 1.10 0.83 1.45
## Zuchwil (Soloth.-Emme) 1.08 0.82 1.41
## Grenchen (Buerenamt) 1.02 0.77 1.36
## Saanen 0.98 0.74 1.30
## Thunersee 0.91 0.70 1.19
## Winznau (Zv Olten) 0.89 0.67 1.18
## Delemont (Sede) 0.88 0.67 1.15
## Fribourg (Aele) 0.87 0.67 1.14
## Moossee-Urtenenbach 0.86 0.65 1.14
## Porrentruy (Sepe) 0.82 0.62 1.07
## Vuippens (Ais) 0.79 0.60 1.05
## Worblental 0.79 0.60 1.05
## Burgdorf 0.78 0.57 1.06
## Aarwangen (Zala) 0.76 0.57 1.00
## Mittl. Emmental 0.60 0.45 0.79
## La-Chaux-de-Fonds 0.58 0.44 0.77
## Region Bern 0.51 0.39 0.67Results are similar in terms of rankings. Of course the reference is now the between-ARA mean instead of one arbitrary ARA, which is more easily interpreted. As expected ARA-level intercepts are pulled towards the mean.
Model A5.3: space-time interaction
We now allow different time trends across ARAs.
if(controls$rerun_models) {
ma5.3 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="iid") +
f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
group=ara2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
data = ww_reg,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.3,file=paste0("../",controls$savepoint,"ma5.3.rds"))
} else {
ma5.3 = readRDS(file=paste0("../",controls$savepoint,"ma5.3.rds"))
}
summary(ma5.3)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 0.798, Running = 110, Post = 0.588, Total = 111
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 13.899 10.975 -9.345 13.899 37.137 13.898 0.005
## weekend -0.001 0.020 -0.040 -0.001 0.037 -0.001 0.000
## hol 0.104 0.058 -0.009 0.104 0.217 0.104 0.000
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 IID model
## day1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.908 0.072 2.770 2.907 3.054 2.903
## Precision for below_loq 0.968 1.195 0.038 0.528 4.969 0.119
## Precision for below_lod 0.005 0.002 0.002 0.005 0.011 0.004
## Precision for day 0.010 0.001 0.009 0.010 0.013 0.010
## Precision for ara1 11.751 2.064 7.697 11.585 15.971 11.784
## Precision for day1 1.166 0.119 0.946 1.161 1.414 1.153
##
## Watanabe-Akaike information criterion (WAIC) ...: 265619.78
## Effective number of parameters .................: 634.31
##
## Marginal log-Likelihood: -159294.52
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')## exp(beta) 0.025quant 0.975quant
## weekend 1.00 0.96 1.04
## hol 1.11 0.99 1.24
## [1] "Relative viral load by ARA compared to average:"ma5.3$summary.random$ara1 %>%
as_tibble() %>%
dplyr::transmute(ara1=ID,
`exp(beta)`=round(exp(mean),2),
`0.025quant`=round(exp(`0.025quant`),2),
`0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>%
dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>%
dplyr::arrange(-`exp(beta)`) %>%
select(-ara1) %>%
column_to_rownames(var="ara_name")## exp(beta) 0.025quant 0.975quant
## Laupen 1.63 1.38 1.92
## Lauterbrunnen 1.52 1.04 2.30
## Grindelwald 1.48 1.25 1.76
## Neuchatel 1.43 1.21 1.70
## Biel 1.24 0.86 1.83
## Interlaken 1.21 0.82 1.80
## Zuchwil (Soloth.-Emme) 1.14 0.96 1.35
## Lyss 1.10 0.75 1.62
## Saanen 1.05 0.73 1.51
## Colombier (La Saunerie) 1.04 0.71 1.52
## Grenchen (Buerenamt) 1.00 0.68 1.46
## Winznau (Zv Olten) 0.96 0.66 1.40
## Thunersee 0.95 0.81 1.13
## Vuippens (Ais) 0.93 0.64 1.34
## Burgdorf 0.92 0.56 1.51
## Moossee-Urtenenbach 0.91 0.62 1.34
## Worblental 0.90 0.62 1.31
## Fribourg (Aele) 0.88 0.75 1.05
## Aarwangen (Zala) 0.87 0.59 1.28
## Porrentruy (Sepe) 0.85 0.72 1.01
## Delemont (Sede) 0.82 0.70 0.97
## La-Chaux-de-Fonds 0.72 0.48 1.05
## Mittl. Emmental 0.70 0.47 1.02
## Region Bern 0.52 0.44 0.62ndays = length(unique(ww_reg$day))
nara = length(unique(ww_reg$ara1))
ma5.3$summary.random$day1 %>%
bind_cols(day=rep(0:(ndays-1),nara),
ara1=rep(1:nara,each=ndays)) %>%
left_join(corr_reg,by = join_by(ara1)) %>%
ggplot() +
geom_hline(yintercept=1,linetype=2,alpha=.5) +
geom_ribbon(aes(x=day,ymin=exp(`0.025quant`),ymax=exp(`0.975quant`),fill=ara_name),alpha=.5) +
geom_line(aes(x=day,y=exp(mean),colour=ara_name)) +
facet_wrap(~ara_name) +
scale_colour_discrete(guide="none") +
scale_fill_discrete(guide="none") +
scale_y_continuous(trans="log",breaks = c(.1,1,10)) +
coord_cartesian(ylim=c(.05,20)) +
labs(title="Deviations from average time trend by ARA",x="Day",y="Relative viral load by ARA") 
This brings a clear improvement in WAIC. We observe similar average between-ARA heterogeneity, with highest relative viral loads in Laupen, Lauterbrunnen, Grindelwald and Neuchâtel (four ARAs where the credible intervals are the highest) and lowest relative viral loads in Delemont, La-Chaux-de-Fonds, Emmental and Bern (four ARAs where the credible intervals are the lowest). While time trends are generally aligned, we observe deviations from the regional average time trend in some ARAs such as Colombier, Delemont and La-Chaux-de-Fonds (comparatively lower viral loads during Summer 2022), and in Lauterbrunnen, Grindelwald and Interlaken (comparatively higher during Summer 2022).
Model A5.4: spatial structure
We now consider the neighboring structure between ARAs.
# setup neighboring matrix
shapes_reg = shapes$ara_shp %>%
dplyr::filter(ara_id %in% ww_reg$ara_id) %>%
left_join(corr_reg,by = join_by(ara_id)) %>%
dplyr::arrange(ara1)
sf_use_s2(FALSE)
graph_reg = spdep::poly2nb(shapes_reg)## although coordinates are longitude/latitude, st_intersects assumes that they are planarpath_graph = paste0("../",controls$savepoint,"W_reg_",select_NUTS2,".adj")
nb2INLA(path_graph, graph_reg)
if(controls$rerun_models) {
ma5.4 = INLA::inla(vl ~ 1 +
f(below_loq,model="iid") +
f(below_lod,model="iid") +
f(day,model="rw2", scale.model=TRUE, constr=TRUE,
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))) +
f(weekend,model="linear",mean.linear=0,prec.linear=.2) +
f(hol,model="linear",mean.linear=0,prec.linear=.2) +
f(ara1,model="bym2",
graph=path_graph,
scale.model = TRUE, constr = TRUE,
hyper = list(theta1 = list("PCprior", c(1, 0.01)), # Pr(sd<1) = 0.01, unlikely to have rr>3just based on the spatial confounding
theta2 = list("PCprior", c(0.5, 0.5))) # Pr(phi<0.5)=0.5, we state that we believe that the unmeasured spatial confounding is driven 50% from the structured and 50% from the unstructured random effect
) +
f(day1,model="rw1", scale.model=TRUE, constr=TRUE,
group=ara2, control.group=list(model="iid"),
hyper=list(prec = list(prior = "pc.prec", param = c(1, 0.01)))),
data = ww_reg,
family = "gamma",
control.compute = list(waic=TRUE,config=TRUE),
control.predictor = list(compute=TRUE,link=1))
saveRDS(ma5.4,file=paste0("../",controls$savepoint,"ma5.4.rds"))
} else {
ma5.4 = readRDS(file=paste0("../",controls$savepoint,"ma5.4.rds"))
}
summary(ma5.4)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", "
## blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 8.61, Running = 142, Post = 1.54, Total = 152
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 13.908 12.669 -11.510 13.908 39.319 13.908 0.002
## weekend -0.002 0.020 -0.040 -0.002 0.037 -0.002 0.000
## hol 0.104 0.058 -0.009 0.104 0.217 0.104 0.000
##
## Random effects:
## Name Model
## below_loq IID model
## below_lod IID model
## day RW2 model
## ara1 BYM2 model
## day1 RW1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision parameter for the Gamma observations 2.913 0.073 2.767 2.914 3.056 2.919
## Precision for below_loq 3.816 3.145 0.474 2.971 12.057 1.372
## Precision for below_lod 0.005 0.008 0.000 0.003 0.025 0.001
## Precision for day 0.011 0.002 0.007 0.010 0.016 0.010
## Precision for ara1 10.756 4.008 5.168 10.005 20.694 8.673
## Phi for ara1 0.183 0.182 0.008 0.118 0.691 0.020
## Precision for day1 1.202 0.148 0.949 1.188 1.533 1.154
##
## Watanabe-Akaike information criterion (WAIC) ...: 265628.62
## Effective number of parameters .................: 638.86
##
## Marginal log-Likelihood: -159271.80
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')## exp(beta) 0.025quant 0.975quant
## weekend 1.00 0.96 1.04
## hol 1.11 0.99 1.24
## [1] "Relative viral load by ARA compared to average:"ma5.4$summary.random$ara1 %>%
as_tibble() %>%
dplyr::transmute(ara1=ID,
`exp(beta)`=round(exp(mean),2),
`0.025quant`=round(exp(`0.025quant`),2),
`0.975quant`=format(round(exp(`0.975quant`),2),scientific=FALSE)) %>%
dplyr::left_join(select(corr_reg,ara1,ara_name),by = join_by(ara1)) %>%
dplyr::arrange(-`exp(beta)`) %>%
dplyr::select(-ara1) %>%
dplyr::filter(!is.na(ara_name)) %>%
column_to_rownames(var="ara_name")## exp(beta) 0.025quant 0.975quant
## Lauterbrunnen 1.64 1.10 2.47
## Laupen 1.62 1.38 1.92
## Grindelwald 1.48 1.25 1.76
## Neuchatel 1.43 1.22 1.69
## Interlaken 1.30 0.87 1.95
## Biel 1.28 0.85 1.94
## Zuchwil (Soloth.-Emme) 1.13 0.96 1.34
## Lyss 1.07 0.73 1.58
## Saanen 1.05 0.72 1.53
## Grenchen (Buerenamt) 1.00 0.68 1.48
## Colombier (La Saunerie) 0.99 0.67 1.46
## Thunersee 0.95 0.81 1.12
## Winznau (Zv Olten) 0.95 0.64 1.41
## Vuippens (Ais) 0.91 0.62 1.34
## Fribourg (Aele) 0.88 0.75 1.04
## Moossee-Urtenenbach 0.86 0.59 1.27
## Burgdorf 0.85 0.51 1.41
## Porrentruy (Sepe) 0.85 0.72 1.00
## Worblental 0.84 0.58 1.24
## Aarwangen (Zala) 0.83 0.56 1.24
## Delemont (Sede) 0.82 0.69 0.97
## La-Chaux-de-Fonds 0.70 0.46 1.03
## Mittl. Emmental 0.67 0.45 0.98
## Region Bern 0.52 0.44 0.62ndays = length(unique(ww_reg$day))
nara = length(unique(ww_reg$ara1))
ma5.4$summary.random$day1 %>%
bind_cols(day=rep(0:(ndays-1),nara),
ara1=rep(1:nara,each=ndays)) %>%
left_join(corr_reg,by = join_by(ara1)) %>%
ggplot() +
geom_hline(yintercept=1,linetype=2,alpha=.5) +
geom_ribbon(aes(x=day,ymin=exp(`0.025quant`),ymax=exp(`0.975quant`),fill=ara_name),alpha=.5) +
geom_line(aes(x=day,y=exp(mean),colour=ara_name)) +
facet_wrap(~ara_name) +
scale_colour_discrete(guide="none") +
scale_fill_discrete(guide="none") +
scale_y_continuous(trans="log",breaks = c(.1,1,10)) +
coord_cartesian(ylim=c(.05,20)) +
labs(title="Deviations from average time trend by ARA",x="Day",y="Relative viral load by ARA") 
Results are very similar to model A5.3. We find that the spatial structure explains between 0 and 70% of the spatial heterogeneity.